Phase transformation in Si from the semiconducting diamond to the metaUic /3-Sn 
phase in QMC and DFT under hydrostatic and anisotropic stress. 
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Silicon undergoes a phase transition from the semiconducting diamond phase to the metallic /3-Sn 
phase under pressure. We use quantum Monte Carlo calculations to predict the transformation pres- 
sure and compare the results to density functional calculations employing the LDA, PBE, PW91, 
WC, AMOS, PBEsol and HSE06 exchange-correlation functionals. Diffusion Monte Carlo predicts a 
transition pressure of 14.0 ± 1.0 GPa slightly above the experimentally observed transition pressure 
range of 11.3 to 12.6 GPa. The HSE06 hybrid functional predicts a transition pressure of 12.4 GPa 
in excellent agreement with experiments. Exchange-correlation functionals using the local-density 
approximation and generalized-gradient approximations result in transition pressures ranging from 
3.5 to 10.0 GPa, well below the experimental values. The transition pressure is sensitive to stress 
anisotropy. Anisotropy in the stress along any of the cubic axes of the diamond phase of silicon 
lowers the equilibrium transition pressure and may explain the discrepancy between the various ex- 
perimental values as well as the small overestimate of the quantum Monte Carlo transition pressure. 

PACS numbers: 64.70.K-, 71.15.Mb, 02.70.Ss 
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I. INTRODUCTION 

Phase transformations between insulating or semicon- 
ducting and metallic phases present a challenge for many 
current theoretical methods. These transitions provide 
a testing ground for comparing the accuracy of quan- 
tum Monte Carlo (QMC) methods and novel exchange- 
correlation functionals used in density functional theory 
(DFT). In this work, we investigate the diamond to /3- 
Sn phase transition in silicon with QMC methods and 
with various types of exchange-correlation functionals in 
DFT. 

Under pressure, Si displays eleven phases with a steady 
increase in coordination number and a transition from 
semiconductor to metal. ^ At ambient pressure. Si occurs 
in the diamond structure. At a pressure of about 12 GPa, 
Si transforms to the high-pressure /3-Sn structure.^ This 
phase transition coincides with a semiconductor-to-metal 
transition and an increase in the coordination number 
from 4 to 6. The transformation reduces the space group 
symmetry from cubic Fd3m (227) to tetragonal I4i/amd 
(141). At the transition pressure, the volume decreases 
by 21%.^'"' It is noteworthy that the melting of Si at am- 
bient conditions also displays a transition both from four- 
fold to sixfold coordination and from the semiconducting 
solid phase to a metallic liquid. 

The observed transition pressure in Si from the dia- 
mond to the /3-Sn phase depends strongly on the ex- 
perimental conditions and is affected by non-hydrostatic 
stresses. Under nominally hydrostatic conditions in dia- 
mond anvil cell experiments using a pressure medium, 
the observed transition pressure ranges from 11.3 to 
12.6 GPa.^"^ Without a pressure medium or under uniax- 
ial compression. Si transforms at significantly lower pres- 



sures of 8-9 GPa.^'^" In shock compression, the transfor- 
mation is observed in a range of 10-14 GPa.^^"^^ 

Previous calculations using DFT determined the 
phase transition pressure for various semilocal exchange- 
correlation functionals. Comparing the results of zero- 
temperature calculations directly with room-temperature 
experiments requires taking into account the phonon con- 
tributions to the free energy. Zero-point energy and finite 
temperature phonon entropy contributions to the free en- 
ergy lower the transition pressure by 1.0 and 0.3 GPa, 
respectively."'^'' Taking this into account, the previous hy- 
drostatic compression calculations gave transition pres- 
sures of 5.7-6.7 GPa using the local density approxima- 
tion (LDA) and 10.1-10.9 GPa using the PW91 func- 
tional, a generalized gradient approximation (GGA) - all 
below the experimental range. '^"''' For non- hydrostatic 
compression, DFT calculations confirm the experimental 
observation that stress anisotropics can lower the transi- 
tion pressure. The transition pressure decreases linearly 
with increasing deviatory stresses along one of the cubic 
axes."'^* 

In this study, we perform QMC and DFT calculations 
to benchmark semilocal and hybrid exchange-correlation 
functionals for the diamond to /3-Sn transformation in 
Si and determine how non-hydrostatic stress affects the 
transition pressure. Section II introduces the DFT and 
QMC methods. We analyze the accuracy of the vari- 
ous approximations in our QMC calculations to deter- 
mine the accuracy of the transition pressure prediction. 
Section III compares the results of the QMC and DFT 
calculations for the transition pressure between the two 
Si phases with the experimental results. The LDA'^ 
predicts a transition pressure that is too low. The 
GGA functionals PBE^" and PW91^^ improve the pre- 
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diction but still underestimate the pressure. The more 
recently developed GGA functional WC^^, AM05^^ and 
PBEsol^'' perform worse than the PBE and PW91 func- 
tionals and predict transition pressures more similar to 
the LDA. The hybrid functional HSEOe'^^'^^ gives a tran- 
sition pressure of 12.4 GPa, in excellent agreement with 
the range of experimental values of 11.3-12.6 GPa. Diffu- 
sion Monte Carlo (DMC) predicts a transition pressure of 
14.0 ±1.0 GPa, which is slightly higher than the HSE06 
results and the range of experimental values. We discuss 
the origin of the shortcomings for the LDA and GGA 
functionals. For the LDA, the lack of gradient terms of 
the functional results in an underestimation of the en- 
ergy difference between the diamond phase and the more 
homogeneous /3-Sn phase. Based on a similar underesti- 
mation of the formation energies for interstitial defects 
in Si,^^'^^ we argue that the lack of nonlocal exchange 
in the LDA and semilocal GGA functionals is responsi- 
ble for the remaining underestimation of the energy dif- 
ference between the semiconducting diamond and /3-Sn 
phases. ^'''^^ 

From the dependence of the energy on volume and c/a 
ratio, we determine the effect of stress anisotropy on the 
transition pressure. We show that stress anisotropics, 
which may be present under experimental loading condi- 
tions, lower the transition pressure and may explain the 
broad range of experimental values. Section IV summa- 
rizes the results. 



II. METHODS 

To accurately compare the predictions of various 
exchange-correlation functionals in DFT to QMC calcu- 
lations, we have to either eliminate or control the physi- 
cal and numerical approximations in both computational 
methods. For the DFT calculations we eliminate or con- 
trol all approximations (other than the choice of the ap- 
proximate exchange-correlation functional) to determine 
the transition pressure accurate to within 0.2 GPa. In 
our QMC calculations we reduce the error introduced by 
the controlled approximations to result in an uncertainty 
of the transition pressure of 1 GPa. In the following, we 
describe our DFT and QMC approaches and estimate the 
accuracy of the approximations. 



A. Density functional methods 

The three main approximations for the density func- 
tional calculations, besides the choice of exchange- 
correlation functional, are the description of the core 
electrons, the accuracy of the basis set and the Bril- 
louin zone integration. For DFT calculations using 
semilocal exchange-correlation functionals, we explicitly 
include the core-electrons in the calculation and avoid 
the pseudopotential approximation by performing all- 
electron calculations using the LAPW method imple- 
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FIG. 1. (color online) Accuracy of the frozen-core PAW 
method. The energy of the diamond and /3-Sn phases of sil- 
icon are shown for the all-electron LAPW method (Wien2k) 
and the frozen core PAW method (VASP) using the PBE 
functional . The energies are shifted such that the minimum 
energy of the diamond phase is at zero. The differences in 
energy predicted by the two methods at any given volume 
are within 5 meV/atom for the diamond phase and within 
10 meV/atom for the /3-Sn phase. The transition pressures 
agree within 0.2 GPa. 

mented in the Wien2k code.^^ The parameters of the 
LAPW basis and the Brillouin zone integration are cho- 
sen to achieve a total energy accuracy of 1 meV/atom. 
This requires a muffin-tin radius Rmt = 1.9 oq and a 
value of Rmt fcmax = 8.0, where fcmax is the planewave 
cutoff of the basis. The Brillouin zone integration is per- 
formed on a 9 X 9 X 9 fc-point mesh for the semiconduct- 
ing diamond phase and al5xl5xl5 /c-point mesh for 
the metallic /3-Sn phase. The calculations for the hy- 
brid functional HSE06 are performed using the VASP 
program (Vienna ab-initio simulation program) employ- 
ing the PAW method within the frozen core approxima- 
tion. ■^"'•^^ A cutoff energy of 500 eV and a 15 x 15 x 15 
/c-point mesh ensure convergence of the total energy to 
1 meV/atom. 

We test the accuracy of the frozen core PAW approxi- 
mation by comparing calculations for the PBE functional 
with the all-electron LAPW method. Figure 1 illustrates 
the close agreement between the frozen core PAW and 
the LAPW method for the diamond and /3-Sn phases of 
silicon. At any given volume, the two methods agree to 
better than 10 meV/atom and the transition pressures 
agree to within 0.2 GPa. 



B. Quantum Monte Carlo methods 

There are two forms of QMC methods that are com- 
monly used for electronic structure calculations, the sim- 
pler variational Monte Carlo (VMC) and the more so- 
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TABLE I. Convergence of controlled approximations in QMC. 
The parameters of the calculations are chosen to reduce the 
errors introduced by the controlled approximations to below 
20 meV/atom. 



Controlled 


Convergence 


Parameter 


approximation 


[meV/atom] 


value 


Statistical error 


<3 


30,000 steps/iVatom 


Orbital interpolation grid 


< 5 


4.5-6.5 points/A 


DMC population control 


<1 


1,000 walkers 


DMC time step 


< 5 


T = 0.025 Ha"^ 


Finite size 


< 20 


up to 128 atoms 



phisticated diffusion Monte Carlo (DMC) method.^^'^^ 
VMC calculates quantum mechanical expectation val- 
ues using Monte Carlo techniques to evaluate the many- 
dimensional integrals. Accuracy of the VMC results de- 
pends crucially on the quality of the trial wave function, 
but DMC can remove most of the error in the trial wave 
function. DMC is a stochastic projector method that 
projects out the ground state from the trial wave func- 
tion. To avoid the fermion sign problem, one typically 
imposes the boundary condition that the nodes of the 
many-body wave function are the same as those of a trial 
wave function. The resulting error, known as the fixed- 
node error, can be reduced by improving the trial wave 
function^'''^^ 

The QMC calculations are performed using the 
CHAMP code.'^^ A norm-conserving Hartree-Fock pseu- 
dopotential eliminates the Is, 2s and 2p electrons of Si 
from the calculation.'^^ The QMC trial wave function con- 
sists of a product of a single Slater determinant of DFT 
orbitals and a Jastrow correlation factor. The orbitals 
of the Slater determinant come from a DFT calculation 
with the CPW2000 code of J.-L. Martins using the PBE 
functional. For the diamond phase the L-point was cho- 
sen to reduce the finite-size error. For the /3-Sn phase the 
(1/2,0,0) fc-point was selected to reduce the finite-size 
error and to avoid any fractional occupancies of the or- 
bitals. The Jastrow factor describes the electron-electron 
and electron-nuclei correlations. Minimizing the energy 
in VMC optimizes the parameters of the Jastrow. •^^''^^ Fi- 
nally, DMC calculations using the optimized trial wave 
function determine the energy of the phases. 

Approximations in the QMC calculations can be classi- 
fied into controlled approximations with systematically- 
reducible error and uncontrolled approximations whose 
errors are unknown and are not systematically reducible. 

Convergence of the controlled approximations. Con- 
trolled approximations include the statistical error of the 
Monte Carlo method, the interpolation grid for the nu- 
merical orbitals, the system size, the number of config- 
urations (walkers), and the time step in the diffusion 
Monte Carlo calculation. Table I summarizes the accu- 
racy of the controlled approximations in the QMC cal- 
culations. All the errors introduced by the controlled 
approximations are reduced such that the resulting tran- 
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FIG. 2. (color online) Finite-size extrapolation of the DMC 
energies for the diamond and /3-Sn phases of Si. The zero of 
the energy is taken to be the extrapolated value for the di- 
amond structure including finite-size corrections. The finite- 
size corrections using the finite-size exchange-correlation func- 
tional method'^® greatly reduce the finite-size error and enable 
an accurate finite-size extrapolation using the corrected DMC 
energies of the 16, 54 and 128 atom cells. 



sition pressure is accurate to within 1 GPa, corresponding 
to an energy accuracy of about 20 meV/atom. 

QMC calculations for extended systems require a 
finite-size extrapolation. The periodic boundary condi- 
tions employed in the calculation lead to artificial corre- 
lations of the electrons. Several methods have been de- 
veloped to reduce the finite-size error such as the model 
periodic Coulomb potential, the completion of the 
structure factor'*^ and the finite-size exchange-correlation 
functional approach. Here we use the last approach, 
which employs a local density functional fit to DMC cal- 
culations for finite system sizes. The finite size cor- 
rections are evaluated using the implementation of the 
finite-size exchange-correlation functional in the PWSCF 
code.^'^ 

Figure 2 shows the DMC energies for the diamond and 
/3-Sn structure of Si at their respective equilibrium vol- 
umes as a function of the number of atoms. Extrapo- 
lation from the raw DMC energies is difficult at best. 
Using the finite-size exchange-correlation functional cor- 
rection greatly improves the extrapolation of the energy 
to infinite system size, particularly for the /3-Sn phase. 
■^^ Using calculations with up to 128 atoms in the unit 
cell results in an extrapolation error below 20 meV/atom 
leading to an uncertainty in the predicted transformation 
pressure of about 1 GPa. 

Uncontrolled approximations. The uncontrolled ap- 
proximations in QMC include the pseudopotential ap- 
proximation for the core electrons, the locality approxi- 
mation for the evaluation of the nonlocal pseudopotential 
terms and the fixed- node approximation in DMC. 

Alfe et al.^'^ showed that core-polarization, which is 
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FIG. 3. (color online) Comparison of the predicted equations 
of state of silicon for diffusion Monte Carlo energies and the 
LDA, PBE and HSE06 exchange-correlation functionals. The 
energy is given relative to the energy of the diamond phase at 
a volume of 20 A^. The DMC energy curve includes the core- 
polarization correction of 30 meV/atom described in Ref. 44. 
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FIG. 4. (color online) Comparison of predictions of the transi- 
tion pressure for various exchange correlation functionals with 
DMC and experiment. The LDA approximation (red) and all 
tested GGA approximations (green), underestimate the tran- 
sition pressure. The HSE06 hybrid functional (blue) agrees 
with the experiment. DMC (black) slightly overestimates the 
transition pressure. The shaded region indicates the the esti- 
mate of the uncertainty of the DMC transition pressure from 
the controlled approximations. 



not included in our QMC calculations, is important for 
the energy difference between the /3-Sn and the diamond 
phase and lowers the energy difference by 30 meV/atom. 
To account for this effect, we apply a constant energy 
shift of 30 meV/atom to the DMC energy of the /3-Sn 
phase. 

In DMC, the trial wave function is used to evaluate the 
non-local part of the pseudopotential. This pseudopo- 
tential locality approximation leads to a non- variational 
error in the DMC energies. Empirically, this error in- 
troduced by the pseudopotential locality approximation 
is usually quite small for well-optimized trial wave func- 
tions.^^ 

The fixed-node error is difficult to estimate and could 
affect the results of our calculation. Possible methods 
to improve the nodes of the trial wave function include 
orbital optimization^^''*'^ and backflow transformation.'^^ 
Both approaches are computationally very demanding 
and beyond the scope of this work. Recent DMC cal- 
culations for self-interstitials in Si'** provide an estimate 
of the size of fixed-node error. These calculations show 
that the backflow transformation significantly improves 
the accuracy of the trial wave function; the variance is 
reduced threefold by the backfiow transformation. The 
total energy for the diamond Si phase is only reduced 
by 12 meV/atom. Furthermore, a partial error cancel- 
lation is observed for differences in total energy between 
the perfect bulk and the defects. We expect a similar 
partial error cancellation of the fixed-node error between 
the diamond and the /3-Sn phase. 



III. RESULTS 

A. Comparison between quantum Monte Carlo 
and various exchange-correlation functionals 

Figure 3 compares the predictions of the different 
exchange-correlation functionals for the energy as a func- 
tion of volume of the diamond and /3-Sn phases of silicon 
with the results of the diffusion Monte Carlo calculations. 
For clarity, we show only the results for the LDA, PBE 
and HSE06 functional. 

Table II presents the equilibrium volume Vq , bulk mod- 
ulus B and pressure dependence of the bulk modulus B' 
for Si in the diamond and /3-Sn phases for the various the- 
oretical methods. For the diamond phase, we fit a Birch- 
Murnaghan equation of state*^ to the energy as function 
of volume over a range from 16.0 to 22.0 A^/atom. Com- 
paring the results with experiment shows an excellent 
agreement of the structural properties with experimental 
values for the more recent GGA functionals, WC, AMOS 
and PBEsol, and the hybrid functional HSE06. The 
LDA, as usual, underestimates the volume, and the GGA 
functionals PBE and PW91 overestimate it. A similar 
agreement is also observed for the volume Vt at the ex- 
perimental transition pressure of 11.7 GPa of Ref. 4. All 
semilocal functionals somewhat underestimate the bulk 
modulus with the largest discrepancy observed for the 
PBE and PW91 functional. All functionals reproduce 
the pressure dependence of the bulk modulus very well. 
Within statistical accuracy, the DMC results agree with 
the experimental values for the diamond phase, especially 
for the equilibrium volume Vq ■ The somewhat large error 
bars on B and B' make a comparison difficult. 
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TABLE II. Equation of state of the diamond and /3-Sn phases of silicon in DMC and DFT with various functionals. Shown 
are the cqxiilibrium volume, Vq, the bulk modulus, B, and the pressure dependence of the bulk modulus, B' for both phases 
and the co/oo ratio for the /3-Sn phase. For comparison with experiment, we list the volume, Vt and the ct/at ratio at the 

experimental transition pressure of 11.7 GPa of Ref. 4. The experimental equilibrium volume of the diamond phase and the 
elastic properties are from Ref. 50. In addition, AEo, denotes the energy difference between the minima of the two phases. 
The transition pressure, pt includes zero-point and finite-temperature corrections which lower the transition pressure by 1.0 
and 0.3 GPa, respectively.^'' The QMC transition pressure also includes a core-polarization correction which lowers the energy 
difference between the diamond and /3-Sn phase by 30 meV/atom.** 
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* At the experimental transition pressure of Ref. 4. 

^ The c/a ratio is not optimized in DMC and fixed at the experimental c/a = 0.550. 
° A value of B' = 4.6 is assumed in the DMC Birch-Murnaghan equation of state. 



For the /3-Sn phase, we calculate the energy as a func- 
tion of volume V and c/a ratio. For each volume, we 
fit a cubic polynomial to determine the minimum en- 
ergy and corresponding c/a ratio. Similar to the analysis 
for the diamond phase, a Birch-Murnaghan equation of 
state"*^ is then fit to the resulting energies as a function 
of volume ranging from 13.2 to 17.1 A^/atom. Since the 
/3-Sn phase transforms under decompression into the R8 
phase, no experimental values are available for the equi- 
librium crystal structure. Instead we compare the volume 
and c/a ratio at the experimental transition pressure of 
11.7 GPa of Ref. 4. The LDA and AM05 functional un- 
derestimate the volume at that pressure while all other 
functionals agree quite well with the experiment. 

The transition pressure pt is determined from the equa- 
tion of state using the common-tangent rule. To com- 
pare the calculated bare transition pressure with exper- 
imental values, we need to include the effect of zero- 
point vibrations and finite temperature, which differ 
in the two phases. -"^^ The results for pt shown in Ta- 
ble II include zero-point and finite-temperature correc- 
tions, which lower the transition pressure by 1.0 and 
0.3 GPa, respectively.-'^^ Figure 4 compares the pre- 
dictions for the transition pressure among the various 
exchange-correlation functionals and DMC with exper- 
iments. Somewhat surprisingly, all semilocal hmction- 
als underestimate the transition pressure, and the GGA 



functionals provide a large range of pressure predictions 
from 3.5 to 10.0 GPa. The PW91 functional predicts 
a pressure of 10.0 GPa close to the experimental range 
from 11.3 to 12.6 GPa. The more recent GGA function- 
als, WC, AMOS and PBEsol, significantly underestimate 
the transition pressure, similarly to the LDA. The hybrid 
functional HSE06, however, predicts a transition pressure 
of 12.4 GPa, in excellent agreement with experiments. 

The DMC calculations predict a transition pressure of 
14.0 ± 1.0 GPa. This value includes the same zero-point 
energy and finite-temperature phonon entropy contribu- 
tions to the free energy^^ as the DFT values, and the 
core-polarization correction by Alfe et al/^. The error 
bar of 1 GPa is an estimate of the combined uncertain- 
ties of the controlled approximations. 

Our DMC transition pressure of 14.0 ±1.0 GPa is 
slightly above the experimental range from 11.3 to 12.6 
GPa. The result is lower by 2.5 GPa than the DMC 
result of Alfe et al.^"^ of 16.5 GPa and agrees within er- 
ror bar with the AFQMC result by Purwanto et al.^^ of 
12.6 ±0.3 GPa. The latter work estimated the transition 
pressure from the energy difference at the experimental 
transition pressure, limiting any estimate of the accuracy 
of the AFQMC method to the pressure value. 

The work by Alfe et al.^'^ and our study determined 
the equation of state of both phases. Both DMC calcula- 
tions accurately predict the structural and elastic proper- 
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ties of the diamond phase demonstrating their accuracy. 
The two DMC studies differ in their transition pressure 
prediction. The higher transition pressure predicted by 
Alfe et al. originates from a larger energy difference of 
475 ± 10 meV/atom'*^ versus 424 ± 20 meV/atom in our 
work (both numbers include the core-polarization cor- 
rection). The difference of 2.5 GPa or 50 meV/atom 
between the two DMC results could be caused by differ- 
ences in the pseudopotential locality error, the finite-size 
extrapolation or the fixed-node error. The main differ- 
ences in our approach are the use of energy minimization 
for optimizing the Jastrow factor parameters, a different 
form of the Jastrow and the finite-size extrapolation. An 
improved optimization of the Jastrow parameters by the 
energy minimization method'^^'"^* factor would reduce the 
error of many of the controlled approximations. It would 
also reduce the uncontrolled pseudopotential locality er- 
ror in DMC. In our work we include a finite-size correc- 
tion^^ and perform a finite-size extrapolation at each vol- 
ume for the two phases using the DMC energies for 16, 54 
and 128 atom cells while Alfe et al. use the DMC energy 
of the 128 atom cells without extrapolation and finite-size 
corrections. The combination of pseudopotential locality 
error and different finite-size extrapolation could explain 
the change in the energy difference of 50 meV/atom. 



B. Discussion of the phase stability 

The differences in the transition pressure predictions 
of the various methods are mostly determined by the 
relative phase stability or energy difference between the 
two phases. The good agreement of the HSE06 hy- 
brid functional with the experimental transition pressure 
indicates that the energy difference should be around 
390 meV/atom. The semilocal functionals all give energy 
differences that are too small, while the DMC energy dif- 
ference appears to be a bit too large. In the following, 
we discuss two arguments for the observed ordering of 
the energy difference between the semiconducting four- 
fold coordinated diamond phase and the metallic sixfold 
coordinated /3-Sn phase for the different functionals. 

We first show that the lack of gradient terms in the 
LDA functional^^'^^ results in an underestimate of the 
energy difference between the two phases. Then, we ar- 
gue that in order to accurately predict the energy differ- 
ence between the two phases, the method also needs to 
predict the band gap accurately. Both LDA and GGA 
functionals fail for the band gap. The inclusion of exact 
exchange in the HSE06 functional recovers the band gap 
by improving the derivative discontinuity of the Kohn- 
Sham potential for integer electron numbers. ^^"^^ 

To understand the trend of the various semi-local 
exchange-correlation functionals we note that a similar 
energy ordering for the various functionals and an agree- 
ment between the DMC and HSE06 results also occurs 
for Si single interstitial defects. The interstitial atom 
and its neighboring atoms have an increased coordina- 
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FIG. 5. (color online) Energy landscape of Si as a function 
of volume and c/a ratio for the HSE06 functional. The two 
minima correspond to the diamond and /3-Sn phases. For the 
chosen unit cell, the c/a ratio of the diamond phase is given 
by \/2 instead of 1. The contour lines are shown for every 
50 meV/atom. 



tion number of five or six. For both the Si intersti- 
tial structures and the /3-Sn phase, the increased coor- 
dination results in a more homogeneous charge density 
compared to the diamond phase. This increased homo- 
geneity explains why the LDA functional underestimates 
the energy difference between the two Si phases and the 
Si interstitial formation energies, respectively. The lack 
of any gradient terms in LDA energetically favors more 
homogeneous charge density distributions.^^ GGA func- 
tionals aim to correct this shortcoming of LDA. However, 
PBE and PW91 violate the gradient expansion of Svend- 
sen and von Barth^^ for slowly varying density systems. 
Both functionals have second-order expansion coefficients 
that are too large. One might therefore expect PBE to 
overestimate the effect of inhomogeneity. However, this 
analysis neglects the effect of higher order contributions. 
We observe that the GGA functionals PBE and PW91 
indeed improve the agreement with DMC, HSE06 and 
experiments but the resulting energy difference between 
the phases is still too small. The more recent GGA func- 
tionals, WC, AM05 and PBEsol, which are designed to 
improve the description of solids, result in energy differ- 
ences between the phases at or even below the LDA value 
and do not work for the transition pressure. The gradient 
corrections alone appear insufficient. 

Recent GW calculations by Rinke et al.^^ showed that 
the failure of LDA and GGA functionals for the inter- 
stitial formation energies in Si is due to underestimation 
of the vertical electron affinities of the interstitial defect 
configurations and related to the band-gap problem. One 
might expect that this underestimation of the band gap 
in LDA and GGA also affects the accuracy of LDA and 
GGA functionals for the semiconductor to metal tran- 
sition in Si. Along a path in configuration space going 
from the diamond to the /3-Sn phase, the band gap closes 
and the density of states at the Fermi level increases. 
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FIG. 6. (color online) Effect of stress anisotropy on phase 
transition pressure. The phase transformation from the dia- 
mond to the /3-Sn phase of Si is reduced by anisotropic stress 
along the [001] direction. The transition pressure reduction is 
similar for the various exchange-correlation functionals. For a 
stress anisotropy A along the [001] direction compared to the 
other two cubic directions, the transition pressure is reduced 
by 2.4 A. 

The bandgap closure occurs too early for the LDA and 
GGA functionals compared to HSE06 and the density 
of states at the Fermi level is higher for the LDA and 
GGA functionals as well. This difference in groundstate 
properties predicted by the LDA and GGA functionals 
correlate with an energy difference along the path that 
is too low for the LDA and GGA functionals. The in- 
clusion of short-ranged Hartree-Fock exchange in the hy- 
brid HSE06 functional is optimized to describe covalent 
bonds and improves the band gap in semiconductors.^^ 
This may be the reason for the good agreement of the 
HSE06 results for the transition pressures in Si with ex- 
periments and for the interstitial formation energies with 
GW calculations^^ and DMC.^^ 



C. Effect of stress anisotropy on the transition 
pressure 

The predicted transition pressure in DMC is at the 
upper end of the experimental values. There are var- 
ious possibilities for why the DMC transition pressure 
prediction is somewhat high. DMC might overestimate 
the energy difference between the metallic /?-Sn phase 
and the semiconducting diamond structure due to the 
fixed-node error of the trial wave function. In order to 
obtain the transition pressure accurate to within 1 GPa, 
one has to determine the energy difference between the 
two phases to within 20 meV/atom. Because of the large 
difference in the electronic structure of the two phases, 
metallic versus semiconducting, it is possible that a suffi- 



ciently accurate cancellation of fixed-node error between 
the two phases does not occur. 

Another possibility is that the measured transition 
pressure could be affected by kinetic effects and stress 
anisotropics in the diamond anvil cell experiments. '^'^'^ 
Kinetic or hysteresis effects typically lead to an increase 
in transition pressure over the thermodynamic equilib- 
rium transition pressure. Stress anisotropy in the sample 
can lower the transition pressure by stabilizing the /3-Sn 
phase over the diamond phase. -'^'''^^ 

We determine the effect of non-hydrostatic stresses on 
the transition pressure following the approach by Cheng 
et al}'' The phase transformation from the diamond to 
the /3-Sn phase occurs when 

F2 - Fi -f I^ < 0, (1) 

where W is the work done by the system during the phase 
transformation and -Fi and F2 are the energies of the di- 
amond and the /3-Sn phase, respectively. In our analysis 
we neglect the effect of stress anisotropy on the zero-point 
energy and finite temperature phonon entropy contribu- 
tions to the free energy and simply include those correc- 
tions as determined by Gaal-Nagy et al}'^ We assume 
a uniform non-hydrostatic compression along the cubic 
and tetragonal axes of the two crystal structures with a 
stress tensor, tr, of the form 

(p-A/3 ° \ 

p-A/3 ' (2) 

P+2A/3 / 

where p is the average applied stress and A is a measure 
of the stress anisotropy between the x or y axes and the 
z-axis. For this loading condition, the work W is 

H2) /-(S) A2) 

W^Px I lylzdlx+Py / Izlxdly+Pz / Ixlydl^, 

J(l) J{1) J{1) 

(3) 

where , py and Pz are the three diagonal entries of the 
stress tensor ct, and Ix, ly and Iz are the lattice constants 
of the crystal structures. The virtual work under non- 
hydrostatic loading is path dependent^* so we calculate 
the work along the shortest path, following Ref. 17. 

To determine the initial and final point of the path 
integral for the virtual work, we calculate the energy 
E{V, c/a) of the Si diamond and /3-Sn phases as a func- 
tion of volume and c/a ratio using the LDA, PBE, PW91 
and IISE06 fimctionals. We leave out the WC, AMOS 
and PBEsol functionals since they significantly underes- 
timate the transition pressure. We use a body-centered 
tetragonal unit cell for the simulations with the three 
lattice vectors (a/2 ao, -\/2 ao, 0), (a/2 oq, — -\/2 oq, 0) and 
(0,0, Co). In these lattice vectors, the c/a ratio of the 
diamond phase is y/2 instead of one. Figure 5 illustrates 
the resulting energy landscape for the HSE06 functional. 
The two minima correspond to the diamond and /3-Sn 
phases of Si. From the energy landscape E{V,c/a), we 
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determine the lattice parameters of the two phases as a 
function of applied stress tensor c using the relations: 

__dE c/a dE 
dV ^ V d{c/ay 

Equations 1 and 3 determine the equilibrium transi- 
tion pressures as a function of stress anisotropy A. Fig- 
ure 6 compares the resulting transition pressures for the 
LDA, PBE and HSE06 functionals with the range of 
experimental values as a function of the applied stress 
anisotropy A. All three functionals predict a similar be- 
havior. The anisotropy in the loading condition A lin- 
early reduces the transition pressure by 2.4 A. The linear 
coefficient is similar for all three functionals and close to 
the values determined previously for the LDA and PW91 
functional. ^'"'^^ Considering that the mechanical strength 
of silicon under uniaxial compression at ambient condi- 
tions is about 7 GPa,^^ it may be reasonable to assume 
that deviatory stresses of the order of 0.5 GPa could 
be present in the diamond-anvil cell experiments with a 
pressure medium and even larger deviatory stresses with- 
out a pressure medium. This would lower the observed 
transition pressure by 1.2 GPa compared to perfect hy- 
drostatic compression. This strong influence of the stress 
anisotropy may explain the range of experimentally ob- 
served transition pressures, particularly the differences 
in the diamond anvil cell experiments with and without 
a pressure medium. The effect of stress anisotropy may 
also explain the difference between the experimental val- 
ues and the DMC prediction. 

IV. CONCLUSION 

We performed calculations of the transition pres- 
sure for the high-pressure phase transformation in Si 
from the semiconducting diamond to the metallic /3- 
Sn phase. Comparisons with experimental values 
benchmark the accuracy of DMC methods and vari- 
ous exchange-correlation functionals in DPT. The hy- 
brid functional HSE06 and the DMC method predict 
similar transition pressures with values of 12.4 and 
14.0 ± 1.0 GPa, respectively, while semilocal LDA and 
GGA functionals predict lower transition pressures rang- 
ing from 3.5 to 10.0 GPa. 



The DMC transition pressure is slightly above the ex- 
perimental range of values of 11.3-12.6 GPa while the 
HSE06 functional agrees with the experiments. The 
DMC energies could be affected by fixed-node error which 
could be reduced using the backflow transformation.^^ 
The LDA prediction of 5.8 GPa is too low. The GGA 
functionals PBE and PW91 improve the prediction but 
still underestimate the pressure. The more recently de- 
veloped GGA functionals, WC, AM05 and PBEsol, per- 
form worse than the PBE and PW91 functionals and pre- 
dict transition pressures more similar to the LDA. Com- 
parison with DMC and GW calculations for point defects 
in Si indicate that the underestimation of the transition 
prcssme may be related to the underestimation of the 
bandgap in Si. 

Calculations for anisotropic loading conditions show 
that the experimental transition pressure could be af- 
fected by stress anisotropy. Stress anisotropy can dra- 
matically reduce the transition pressure for the diamond 
to /9-Sn transformation. An anisotropy of only 0.5 GPa 
along any of the cubic axes reduces the transition pres- 
sure by 1.2 GPa, explaining the range of transition pres- 
sures observed in the diamond anvil cell experiments with 
and without a pressure medium and possibly also the dif- 
ference between the experiments and the DMC. 
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